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We study two planar square lattice Heisenberg models with explicit dimerization or quadrumerization of 
the couplings in the form of ladder and plaquette arrangements. We investigate the quantum critical points of 
those models by means of (stochastic series expansion) quantum Monte Carlo simulations as a function of the 
coupling ratio a = J'/ J. The critical point of the order-disorder quantum phase transition in the ladder model is 
determined as a c = 1.9096(2) improving on previous studies. For the plaquette model, we obtain a c = 1.8230(2) 
establishing a first benchmark for this model from quantum Monte Carlo simulations. Based on those values, 
we give further convincing evidence that the models are in the three-dimensional (3D) classical Heisenberg 
universality class. The results of this contribution shall be useful as references for future investigations on 
planar Heisenberg models such as concerning the influence of non-magnetic impurities at the quantum critical 
point. 
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I. INTRODUCTION 

The study of quantum effects in magnetism is an ongoing 
and fascinating part of physics research^ Within this area, 
the low-dimensional 5 = 1/2 Heisenberg antiferromagnet 
plays an eminent role. This is partly because it correctly de- 
scribes aspects of cuprate superconductors and is thus imple- 
mented in nature. Second, the Heisenberg model and varia- 
tions have seen a lot of investigations as toy models where 
quantum fluctuations lead to unexpectedly rich and exotic 
ground-states (such as valence bond solids and valence bond 
liquids). Recent experiments in optical lattices 3 further pro- 
vide the perspective to directly implement those models in a 
pure environment thereby enabling a direct experimental ac- 
cess and comparison between theory and measurements. 

In two-dimensional (2D) Heisenberg models, the Mermin- 
Wagner theorem forbids phase transitions to occur at T + 0, 
yet quantum fluctuations may lead to a transition between 
ground states, for example from an ordered Neel to a disor- 
dered state at zero temperature. Such transitions are termed 
quantum phase transitions. 4,5 One way in which quantum fluc- 
tuations can destroy order is for example provided via frustra- 
tion of bonds (next-nearest-neighbor couplings) or the inclu- 
sion of 4-site interactions^ 

In a second mechanism, competition between locally vary- 
ing nearest-neighbor bonds of the same kind has been iden- 
tified to cause quantum phase transitions, for example by fa- 
voring the formation of spin singlets. An important class of 
models in which the latter mechanism is at work are the so 
called dimerized Heisenberg models (where we also use the 
term for extended models with quadrumerization, etc.), where 
the competition among couplings is explicitly introduced in a 
geometric manner. Apart from their relevance as simple mod- 
els for quantum phase transitions such systems have been in 
recent focus in connection with Bose-Einstein condensation 
of magnons, 7 A prominent example of dimerized models is 
the 5=1/2 bilayer Heisenberg syste m 8 ' 9 ' 10 ' 11 ' 12 which con- 



sists of two LxL layers, where the inter-layer coupling J x can 
be different from the intra-layer coupling J (both couplings 
antiferromagnetic). Competition between J x and / can drive 
a quantum phase transition. 

Due to progress and availability of unbiased and efficient 
methodological scheme*^ some numerical contributions in 
the literature were lately pushing results on those bilayer sys- 
tems to unprecedented accuracy for quantum models, allow- 
ing for very detailed studies in the quantum critical regime. 
Following the high precision study on two bilayer systems by 
Wang et al.^ Hoglund and Sandvik could for example re- 
port on anomalous response 15 of non-magnetic impurities, for 
which an accurate knowledge of the quantum critical point 
was a prerequisite. The overall interest on such impurity based 
questions is growing ) 16 ' 17 ' 18 therefore asking for the general 
availability of more detailed data also in other systems. 

While the level of accuracy has reached a very high quality 
for bilayer systems, this is not equally the case for planar ge- 
ometries. After the seminal simulation of the CaVO lattice by 
Troyer et al. only the coupled ladder model was considered 
in more detail using quantum Monte Carlo (QMC) studies. 
A main result of these investigations was the confirmation of 
the critical exponents predicted by field theory Z±22k 

In an effort to systematically improve and extend these re- 
sults to other planar Heisenberg models, we have recently 
started with a contribution 23 reporting on peculiar and non- 
universal features of a particular dimerized model, called the 
J - J' or staggered mode h 24 ' 25 ' 26 ' 27 i 28 In Ref. our presen- 
tation is based on a detailed scaling analysis at criticality and 
a comparison between several dimerized models including bi- 
layer and planar geometries. As a prerequisite to this compari- 
son, we have also presented new but preliminary results on the 
ladder and plaquette Heisenberg model without showing any 
details of our numerical data nor its data analysis. An in-depth 
study of these models on its own is, however, useful for sev- 
eral reasons. Apart from the aforementioned motivation con- 
cerning impurities, new benchmark results shall be useful for 
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Figure l : (Color online) (a) Visualization of the ladder model on the 
two-dimensional square lattice. The quantum spin (5 = 1 /2) degrees 
of freedom live on a square lattice with different nearest-neighbour 
couplings J and ]' (thin and thick bonds). The lattice bonds corre- 
sponding to couplings /' are denoted as (;', j)' in the Hamiltonian. 
(b) Similar for the plaquette model, favoring quadrumer formation. 
Both systems are studied using periodic boundary conditions. 

thermodynamic considerations in the quantum critical regime 
and for further developing and testing novel algorithms and 
numerical techniques. 

In order to close this existing gap, we consider in this paper 
the critical points of the ladder and plaquette models defined 
by the Hamiltonian 

•H = /^S,S J + /'^S,S / . (1) 

(u> ajy 

Here, S,- = (1/2) (o~ x , cr y , cr z ) denotes the usual spin-1 /2 oper- 
ator at lattice site i, and J and J' the antiferromagnetic cou- 
pling constants defined on bonds <;, j) and (i, j)', respectively. 
The arrangements of the bonds on the square lattice of size 
L in both directions can be seen in Fig. Q] Let us define the 
quantity a = J' /J as the ratio of the two competing couplings. 
For a > a c > 1 the systems will be disordered and gaped due 
to formation of spin-singlets. For a; < a < a c the systems 
possess Neel order and there is no gap. Here a; is some lower 
boundary at which a second transition can take place. In this 
regard, it is interesting to note that long range Neel order even 
for a = 1 was only recently established rigorously. 29 With a c 
we denote the quantum critical point. Throughout this work 
we fix J = 1 and study the transition from the Neel to a disor- 
dered state, when a (or J') is increased^ 7 - 

Let us first summarize some previous work on the sub- 



ject. Early contributions on the ladder model were done by 
Singh et al.r^ who used series expansions to access the crit- 
ical point. Numerically oriented work followed from Katoh 
and Imada^S and was later improved by Matsumoto et al*& in 
a detailed QMC study which had its major objective in study- 
ing the 5=1 case. For S = 1/2, to our knowledge the best 
known value for the critical coupling is taken from that paper 
as a c = 1.9107(2) (which is the inverse of 0.52337(3)), to- 
gether with an estimate of the critical exponent v = 0.71(3). 
The latter result is often used/quoted in favor of 0(3) univer- 
sality based on field theory. The 5 = 1/2 ladder model has 
been further investigated in three dimensions (3D) in con- 
nection with field induced phenomena and Bose condensa- 
tion of magnonsi^ 2 - The effects of random site dilution in 
the dimerized phase were also studied^ Quite generally, the 
coupled ladder model is nowadays often used as a paradig- 
matic model in discussions of quantum phase transitions and 
quantum magnetism* 2 ^ 

Less is known about the plaquette model, which 
was studied before mainly analytically or with series 
expansions! 35,36,37 A recent study on the quadrumerized 
Shastry-Sutherland-model 38 , using mainly exact diagonaliza- 
tion methods, also contains a (hidden) QMC estimate of the 
critical coupling a c « 1.82 for the pure plaquette model. Ad- 
ditionally, the plaquettized model returned into focus using a 
numerical scheme called contractor renormalization (CORE) 
method. 39 Still, it lacks a detailed quantum Monte Carlo in- 
vestigation as presented in this paper. 

The reason to reconsider the ladder model is threefold. 
First, we like to test our algorithm and approach on known 
models. Our second motivation is to complement the descrip- 
tion of the phase transition in the ladder model beyond to what 
was done earlier. This includes the extension to different criti- 
cal quantities, inclusion of corrections in the finite-size scaling 
analysis and calculation of critical exponents not considered 
before. Our aim is also to make the value of v more accurate 
for definite interpretation in favor of 0(3) universality. Lastly, 
a major obj ective is to derive results which we partly presented 
in Ref. 23| as the dimerized ladder model is so similar to the 
staggered model. 

We organize our paper as follows. In Sec. [TT] we shortly 
present our implementation of the QMC method and data- 
analysis approaches. Standard observables used to detect the 
critical point are defined and discussed. A detailed presenta- 
tion of our numerical data with a focus on the critical point 
is given in Sec. [HI] Section [TVl contains a finite-size scaling 
analysis of the critical exponents and a summary is given in 
SecH 



II. SIMULATION METHODS AND FINITE- SIZE 
SCALING 

A. Quantum Monte Carlo simulations 

In this work, we report on simulations based on our imple- 
mentation of the stochastic series expansion (SSE) method by 
Sandvik and Kurkijarvi42 Due to its discrete nature, this QMC 



3 



2.28 r 
2.26 
2.24 - 
2.22 - 
2.20 
2.18 - 
2.16 
2.14 - 
2.12 - 
2.10 - 



10 



100 



Figure 2: (Color online) Convergence test for the plaquette model 
at a system size L = 32 and coupling a = 1.82 displaying the two 
quantities Q2 (upper curve) and (\m~\) (lower curve). Ground-state 
properties are sampled for /? > 100. The staggered magnetization 
was multiplied by 25 for convenience. 



scheme is a convenient and powerful method to implement. 
The central idea of SSE is to sample the series expansion of 
the partition function 



Z = tr[exp(-/W)]=^^ 



(-j3)"(a\<H n \a) 



(2) 



with n being the expansion order, \a) a basis state of the spin 
space, and /3 the inverse temperature. While the original al- 
gorithm used local Metropolis-type updates, major improve- 
ments were achieved by introducing cluster or operator loop 
updates, 41 Our own implementation is based on the directed 
loop 4 * generalization together with additional ideas described 
by Alet et al. .— The recent incorporation of the Wang-Landau 
method 44 into the SSE scheme 4 - allows the use of multihis- 
togram techniques on QMC dat a 46,47 which is useful to ob- 
tain unbiased continuous curves through data points, a feature 
which we use partly in our data analysis. 

In order to access zero temperature properties of the spin 
system, all simulations must be performed at sufficiently large 
/3 so that quantities of interest assume their ground-state value. 
In this contribution this is done in a two stage procedure. For 
a chosen lattice size, we check explicit convergence of ob- 
servables by a f3 doubling approach, i.e., we double y8 until 
quantities agree within error bars. Once a suitable /3 is fixed 
for the chosen size, standard aspect ratio scaling is employed. 
Hence, we fix /3 at lattice size L according to /3l = sL, with s 
being the scale determined in the doubling schemed FigurefJ] 
shows a particular convergence test for a medium sized lat- 
tice (L = 32) indicating ground-state convergence for/? > 100 
for two exemplary observables defined below. This concrete 
test was performed close to the critical point for the plaquette 
model using 4 x 10 5 sweeps. The inverse temperatures used in 
this study are therefore rather large compared to some earlier 
studies. 



B. Observables 

In order to determine the quantum critical point, we look at 
well-known observables. Next to trivial quantities such as the 
average energy per site e, we consider the staggered magneti- 
zation defined by 



(3) 



where the sum runs over all N 
the usual Binder parameters 



L lattice sites, together with 



Q2 



<K) 2 ) 
(Kl> 2 ' 
(K) 4 ) 
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(4) 
(5) 



These parameters are dimensionless and they possess the 
property to cross at the quantum critical point. Note that 
the staggered magnetization and the Binder parameters can be 
determined quite efficiently by averaging over spin represen- 
tations in the operator direction 42 of the SSE representation. 
The brackets (...) therefore signify (m|) = ((ml) op ) con f. 

Second, we study the correlation length £ of the system. We 
employ the standard second-moment approach, which uses 
the structure factors S (q) defined by 



(6) 



with q being a wave vector in Fourier space and r, the vector 
pointing to site ; on the real space lattice. This quantity can be 
efficiently obtained for arbitrary q during the diagonal update, 



as 



tn-\ 

^m q [p]m q [p]* 
V p 



(7) 



where the index p is running over the operator sequence hav- 
ing n non-unit operators. The quantities m q [p] are defined 
at SSE operator slice p as m Cj [p] = 2f S*[p](cos(qrj) - 
i sin(qr,)) and m q [p\* denotes its complex conjugate. The cor- 
relation length is then estimated by 



6-£ 



(8) 



For the anisotropic ladder model we expect g x + g y on the 
square lattice. We found it most useful to look at the correla- 
tion length in y-direction of the system. This choice is arbi- 
trary but somehow motivated from Ref.|23] because % y showed 
good scaling for the staggered Heisenberg model. From stan- 
dard finite-size scaling theory we expect the quantity £ y /L to 
cross for different lattice sizes at the quantum critical point. In 
case of the symmetric plaquette model, an improved estimate 
for the spatial correlation length can be obtained by taking 



(9) 



4 



Lastly, we consider the spin stiffness p s given by 48 
3 

Ps = — (wl + wh , 
4B - 

with w x ,w y being winding numbers defined by 
w A = (N$-NJ)/L A (A = x,y). 



(10) 



(11) 



The symbols Nt and Nj represent the number of operators of 
type SfS~j and SfSj along the /l-direction in the SSE con- 
figuration. The spin stiffness measures the response in free 
energy upon a boundary twist on the staggered magnetization 
(the order parameter field 0) and is also called superfluid den- 
sity in other contexts. At a quantum critical point in a 2D 
system it is expected to scale as p s ~ L d ~ 2 ~ z , where z is the 
dynamical critical exponentJi 4 ^ 



C. Finite-size scaling 

In this paper, we employ a variety of finite-size scaling 
methods to determine various critical quantities from the 
quantum critical point to the critical exponents. To this end, 
we make use of the standard scaling ansatz in the vicinity of 
the critical point 



L (t) = L xlv go{tL" v ) , 



■1M 



(12) 



where v is the critical exponent of the correlation length, A the 
critical exponent of the quantity O, go( x ) the scaling function, 
t the reduced coupling defined by t — (a - a c )/a c , and L the 
lattice size. 

Analysis to ( fT2l was performed in the previous QMC study 
on the ladder model in Ref. |2(J Here, we would like to go 
one step further and take leading corrections to scaling into 
account. Apart from higher order terms 0(1 /L 2 ), the renor- 
malization group (RG) then predicts a scaling of the form 

L (t) = L A ' V [ 80 (tL l/v ) + L-Ogjfl}'*)] , (13) 

where at is the leading correction exponent and ga>(x) another 
scaling function. Writing gu>{x) = c(x)go(x), this becomes 



L (t)=L A/v (l+c(x)L-ngo(x), 



(14) 



with x = tL l l y and a coefficient c(x) depending on x. To zeroth 
order, and for x small we may set c(x) « c = const and arrive 
at the usually employed form 



L (i) = L xlv (\+cL~")g {x). 



(15) 



We consider (fT~5T > as our primary ansatz in the data analysis. 
Note, that in the literature, another ansatz in form of 



L (t) = L xlv {\ + cL-° , )g (tL l/v + dL-f jv ) , 



(16) 



has been discussed which represents an effective approxima- 
tion to ([TBI in the vicinity of the quantum critical pointJiS 
Here, u and <p represent effective corrections, approximating 
the correct RG behavior. In Ref. [IH which is closely related 



to the present paper, the authors employed (fTol i and obtained 
results in excellent agreement with the expectations. Here, 
we will primarily employ ( fTBT l and in some instances compare 
our result to ( fTBT l. In any case, we use this procedure mainly 
to obtain the critical coupling a c *2. We emphasize that final 
results of critical exponents will be given as obtained from or- 
dinary scaling methods at the critical point (x - 0), which are 
described in Sec. [IV] 

Data analysis according to Eq. <TT3T > is known as "data col- 
lapsing". In practice, this can often be achieved by Taylor 
expanding the scaling function go(x) for x — » into a polyno- 
mial of the form 

go(x) = go + g\x + gix 2 + . . . . (17) 

Using this ansatz, relation ( fT2l is turned into 

L (t) = L A ' V (g Q + L xlv gl t + L 2 ' v g 2 t 2 + ...), (18) 

where all free parameters can then be determined by a 
nonlinear-fit of the measured data. The generalization to dT5b 
is obvious. 

We have recently implemented a related method, which 
does not need to make use of Taylor expanding the function 
g (x)£L Using multihistogram techniques, it is possible to di- 
rectly perform a collapse of the data by minimizing the weight 
function 



dx 



2 L (x) - L (x) 



(19) 



where L (x) = L (t)/(L A/v (l + cL~")) and x = tL 1/v . With 

Ol(x) = YjlOl(x)/hl, we denote the average over ni lattice 
sizes. For the quantities Q\, Q2, % y /L, g/L, and p s L we have 
A/v = 0. 

III. SIMULATION RESULTS AND THE CRITICAL POINT 

We performed various simulations on the ladder and pla- 
quette model for lattice sizes specified in Table [jj employing 
the methods described in the last section. All runs where done 
using periodic boundary conditions. The sample size of mea- 
sured data is of the order of 4 x 10 5 for the plaquette model 
and 8 x 10 5 in case of the ladder model, giving an indication 
that the ladder model is somewhat harder to simulate. We typ- 
ically performed 1 x 10 4 sweeps for equilibration. Measure- 
ments were taken every sweep and each sweep constructed 
as many loops as necessary in order to visit 2n vertices in 
the SSE operator expansion on average. A summary of the 
raw data obtained from the simulations is displayed in Fig. [3] 



Table I: Summary of lattice sizes L studied in the simulations. 



model 



lattice sizes L 



ladder 8, 10, 12, 14, 16, 20, 24, 28, 32, 36, 40, 52, 64 

plaquette 8, 10,12, 16,20,24,28,32,36,40,44,48,56,72 
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Figure 3: Behavior of different observables close to the quantum critical points. Different curves correspond to different system sizes (see 
TableUl, where larger slope means larger system size. The left column displays (a) the spin stiffness multiplied by the system size L, (b) the 
correlation length g y divided by system size L, and (c) the Binder parameter Q2 for the ladder model. The plots (d)-(f) in the right column 
show the same quantities obtained for the plaquette model. 



where we show the spin stiffness p s , the correlation length £ 
and the Binder parameter Q2. The left and right panels in 
Fig. [3] distinguish results for the ladder and plaquette model, 
respectively. Evidently, all quantities cross close to an appar- 
ent quantum critical point justifying the scaling assumptions 
for the observables described above. However, clear finite- 
size corrections can be observed for both cases as the crossing 
points for small lattice sizes show large displacements. This 
behavior is expected and in accordance to the data published 
in Ref. 14 Our hope is that those corrections can be described 
by the correction terms included in the scaling ansatz (fl3l l (or 
(fT6]l). Using the raw data, we will now try to extract a precise 
estimate of the quantum critical point. To reach this aim, we 
will follow a two-stage process, starting with an analysis of 
the crossing points followed by a finite-size scaling investiga- 
tion using the collapsing technique. 

This will in principle also give us estimates of critical expo- 
nents but we leave this issue for a more detailed investigation 



in Sec.[IV]using ordinary and well-established methods. 



A. Estimation of the critical point from curve crossings 

Finite-size scaling analysis with scaling functions involv- 
ing many free parameters is a tedious and difficult task due 
to well-known problems of multidimensional nonlinear min- 
imization. Before we attempt to perform a full finite-size 
scaling study using Eq. ( fT3b , we would therefore like to set 
bounds on the possible values of the critical coupling a c . To 
this end, a convenient approach consists in looking at the scal- 
ing of crossing points of curves at L\ and L2 (where L2 = 2 L\) 
for different values of L\. The crossing points are easily ob- 
tained using either the multihistogram method or fitting data 
at L\ and L2 to the simple scaling ansatz in Eq. ( fT2b (us- 
ing polynomial interpolation). Performing this procedure on 
the various observables of Fig.[3](and Q\) yields the plots of 
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Figure 4: (Color online) Crossing points from data at L and 2L of the 
four quantities p s L, g/L, Qi, and Q 2 versus the inverse lattice size. At 
L = oo all curves should meet and define the quantum critical point 
a c . (a) Analysis for the ladder model, (b) Same for the plaquette 
model. 



Fig. [4] which show convergence of the intersection points to 
the quantum critical point in the thermodynamic limit. The 
plots are presented with an x-axis as 1/L, since we do not 
know the correct scaling a-priori. The qualitative behavior of 
the differe nt q uantities toward the critical point is rather simi- 
lar to Ref.[lJ. We find, that for both models the spin stiffness 
has the least finite-size corrections, followed by the correla- 
tion length, and that the normal Binder parameter Q2 shows 
large deviations at small lattice sizes. This is not necessar- 
ily a disadvantage since this often leads to better controlled 
fits. Before performing some fits, however, let us emphasize 
that in case of the staggered model considered in Ref.|23l the 
spin stiffness displayed a qualitatively different convergence 
toward the infinite-volume limit because there the correlation 
length t; y showed less finite-size effects. This proves that p s is 
not always the best quantity. 

Since all quantities give a rather consistent picture in their 
scaling properties we can safely bracket the critical couplings 
from the crossings using the largest available (L, 2L) pair. This 
yields a c e [1.9070, 1.9105] and a c e [1.821, 1.834] for the 
ladder and plaquette cases, respectively. It is tempting to ob- 
tain a more precise estimate from fitting the crossing points to 



Table II: Estimates for the critical point derived from Eq. ( 1201 ) for the 
ladder (top group) and plaquette model (bottom group). 
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2.6(1) 


0.72 



an ansatz due to Binder— 

ac (L,2L) = a c + -^, (20) 

which states that the crossings should normally converge 
faster than L and would indeed show no L-dependence 
at all if a) — 00, i.e. no correction. In this ansatz b is a con- 
stant and we neglected subleading corrections from the "shift" 
term (p. This term can in principle be included^ leading to 
fits which are more difficult to perform. The smooth curves in 
Fig. |4ja) for the ladder model correspond to fits for the cor- 
relation length, the spin stiffness, and the Binder parameters 
Qi and Q 2 , which yield a c = 1.9097(3) %\ a c = 1.9092(6) 
(p s ), or c = 1.9095(5) (Qi), and a c = 1.9093(3) (Q 2 ). They are 
all in agreement within error bars. For the plaquette model 
(Fig. Eb)) we obtain in the same order a c = 1.8232(3), 
a c = 1.8228(4), a c = 1.8238(8), and a c = 1.8229(4), re- 
spectively. All fit results are summarized in Table [II] where 
we additionally give the fitted quantity 1/v + u> and the qual- 
ity of the fits through the chi-squared per degree of freedom 
(X 1 /d.o.f .). Under the assumption that the correlation length 
exponent v « 0.7, we deduce that a> lies roughly in the in- 
terval [0.8, 1.2] for the correlation length and the Binder pa- 
rameters. For the spin stiffness, interestingly, u> seems to be 
smaller. The stiffness thus appears to cross close to the quan- 
tum critical point but has slow convergence towards it. On the 
other hand, the spin stiffness could not be well described by 
Eq. fl20"| i. A similar effect will, in fact, be seen in the analysis 
of Sec. HUB] 

We feel that the critical points obtained above give a fair 
estimate as they agree within error bars. A posteriori, this 
justifies the neglection of (p. Finally, it should be clear, that by 
the same approach other estimates, like v and go, can and have 
been bracketed aiding in the collapse analysis now to come. 

B. The critical point from data collapses 

In the previous section first estimates of the critical points 
were obtained. Next, our goal is to cross-check and possibly 
improve the accuracy by analyzing the data for the full set of a 
values around the crossing points including all lattice sizes in 
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Table U We will therefore now elaborate on the data collapse 
procedure to the scaling ansatz of Eq. ( fT3T >. knowing that we 
have to include subleading corrections terms. In this process 
we will leave all parameters free, since we want to avoid pre- 
occupation about the universality class. Of course we keep in 
mind the bracketing of some important quantities in the pre- 
vious section. Fitting is done using Eq. (T% or ( TT9b . The two 
approaches have been compared and we could not detect a 
noticeable difference in the outcome. We hence use the less 
time consuming approach according to Eq. (T7\ for which a 
fourth-order polynomial for go{x) is employed. 

Due to potential problems with multidimensional fitting, 
the analysis is repeated for at least two different scenarios. 
In a first case we ignore the subleading shift correction, i.e. 
we set (f> — oo (or d — 0) to obtain a first idea of the criti- 
cal coupling, the correlation length exponent v and other pa- 
rameters. We will see that apart from a few exceptions, this 
approach actually describes our data well enough. Next we 
repeat the collapse taking into account possible shift correc- 
tions, described by a finite tp. All fits are repeated multiple 
times including random noise on the starting parameters as 
well as on the raw data. In the latter case, the noise is taken 
to be normal distributed and within the Jackknife- 3 errors cr 
of the original data points. We typically perform 1000 fits for 
each observable. All quoted error bars are then understood 
as being the error bars from this bootstrap 53 procedure. Fig- 
ure |5j a) outlines this procedure and shows that the collapse is 
well behaved. Random starting values converge to a narrow 
collapse region. Figures |3b,c) display histograms of the final 
critical couplings obtained from the bootstrap procedure for 
the ladder and plaquette model, respectively. It is seen that 




Figure 5: (Color online) (a) Visualization of the collapsing analy- 
sis. The (red) crosses (x) signify the starting values and dense (blue) 
points (+) the final values of (a c , v) of the procedure for the spin 
stiffness in case of the plaquette model. Histograms of the final value 
for ff c for several observables for (b) the ladder model, and (c) the 
plaquette model. 



Table III: Tabulated results for the critical coupling ratio a c , the ex- 
ponent v, and the factor go from the collapse procedure for both the 
ladder (upper group) and the plaquette model (lower group). In some 
cases results from two fits, with and without a <p term are given. 





restr. 


a c 


V 


go 


Qi 


d = 





1.9094(3) 


0.717(10) 2.32(1) 


Q\ 


d = 





1.9096(4) 


0.72(1) 


1.451(8) 




no 




1.9094(3) 


0.72(1) 


1.449(3) 


p s L 


no 




1.90974(15) 0.705(7) 


1.155(10) 


ty/L 


d = 





1.9098(4) 


0.715(10) 0.62(1) 


Qi 


d = 





1.8228(4) 


0.716(6) 


2.313(6) 




no 




1.8227(4) 


0.72(1) 


2.311(5) 




d = 





1.8238(6) 


0.72(1) 


1.453(2) 




no 




1.8228(6) 


0.72(1) 


1.447(4) 




d = 





1.8230(3) 


0.67(1) 


1.28(3) 




no 




1.8230(2) 


0.707(6) 


1.27(2) 




d = 





1.8232(2) 


0.709(6) 


0.706(5) 




no 




1.8231(2) 


0.713(6) 


0.70(1) 



the results are consistent as they more or less overlap, yet we 
note a systematic effect as the Binder parameter tends to give 
smaller estimates in comparison to the correlation length and 
the spin stiffness. This is also in accordance with the data on 
the full bilayer of Ref. [l4|. Table ITII1 summarizes concrete re- 
sults for the different models and observables. The best results 
for a c are obtained from the spin stiffness which usually inter- 
polates between values from the correlation length and Qo. 

Second, we could not detect a noticeable difference in the 
results if we include a (f> degree of freedom. An exception to 
this observation is the spin stiffness, which showed controlled 
fits only in presence of (which probably acts as a kind of 
stabilizer). This fact agrees with the observation made during 
the analysis of the crossing points above but is presently not 
well understood. The results for the exponent v are consistent 
with 0(3) universality. Finally, typical values for u> are in the 
range of u> 6 [0.8, 1.3], consistent with the previous section. 
In case of the spin stiffness, we obtain wxl.4 and <p/v « 2.5. 
Using these results, concrete data collapses of the original data 
are given in Fig. |6]which display a very good collapse quality. 

In principle, one would need to perform additional inves- 
tigations on the influence of size of the collapsed regime x 
(see Eq. $1%). We have done that partly, but do not attempt 
a detailed extrapolation as we will extract the actual critical 
exponents by a different method. In any case, we believe that 
our estimates for a c are correct beyond doubt as they are con- 
sistently obtained from three independent methods (crossing 
analysis, collapse to ( fT3] l, and collapse to (fT6l)). This also jus- 
tifies the use of the approximations which are present in the 
finite-size scaling ansatz. 

We now state the main result of this section in giving our fi- 
nal estimates for the critical couplings. Since no details about 
systematic errors (e.g. from undescribed correction effects 
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0.90 
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Figure 6: (Color online) Data collapses for the ladder and plaquette 
(wider range) model displaying (a) the Binder ratio Q 2 (where the 
collapse for the ladder model was shifted upwards by 0.05 for better 
visibility), (b) the correlation length ij y respectively and (c) the 
spin stiffness p s . Apart from the special case of the spin stiffness, 
collapses are shown without the the <f> correction term. 



etc.) are known, a plain average of the critical coupling es- 
timates from Qi, p s , and % y is probably the best choice (and 
is the same as a weighted average). This way, our final esti- 
mate is a c = 1.9096(2) and a c = 1.8230(2) for the ladder and 
plaquette models, respectively. In case of the ladder model, 
this result is in slight disagreement with the previous value of 
1.9107(2) in Ref. 

Before we go on, it is interesting to observe from the quanti- 
ties go listed in TablelHTl that both the Binder parameters at the 
crossing point seem to be consistent within error bars among 
the two models, whereas the spin stiffness and the correla- 
tion length clearly do not possess this property but the reader 
should keep in mind that ij y and f are slightly different quan- 
tities. 



IV. SCALING AT CRITICALITY 

Having determined estimates for the critical couplings, we 
now turn to an investigation of the critical exponents. To this 
end, we make use of standard methods of Monte Carlo data 
analysis. Our reason to decouple this investigation from the 
collapse analysis is to get independent and unbiased estimates. 
A fit at a predetermined critical point, secondly, has less de- 
grees of freedom and is hence easier to control. 

Analysis of the exponents is performed using standard re- 
lations and definitions. An established method to obtain the 
correlation length exponent v, is via the slope sq 2 = dQil&a 
of the Binder parameter evaluated at the critical point. Using 
Eq. ( fT3l > we arrive at 



l/v 



(21) 



Other exponents, in particular, fi and rj are calculated from the 
order parameter and the structure factor at criticality as 



<KI> 



-jS/v 



S(n,n)~L 



(22) 



where we assume Lorentz invariance, i.e., z - 1 in S (n, n) ~ 

In order to use Eqs. ( 12Tb . ( l22l we need data at the quantum 
critical point. This can in principle be achieved by perform- 
ing new simulations. Since we have rather good data in the 
vicinity of a c already, we instead choose to compute sg,, m z s , 
and S (n, n) from polynomial interpolation or multihistogram 
reweighting as employed in the last section. We have checked 
the consistency of the two approaches and use the first method 
from here on. Again, a bootstrap with 1000 samples is per- 
formed on top of this interpolation, varying the raw input data 
within the uncertainties. Figure [7] summarizes and displays 
the critical data so obtained. All plots are in a In - In style 
vs. the lattice size L. It is evident that straight lines represent 
the data rather well. To make this statement more quantitative 
we now perform and present detailed fits and their results in 
Table |IV| For each quantity, 3 fits are performed correspond- 
ing to the best estimate of a c , as well as its lower and upper 
bounds from the uncertainty. In case of the ladder model, we 
also try a further fit at the previous estimate of Ref. |2(J Sev- 
eral observations can be made regarding our results. First, the 
exponent v obtained from the slope of the Binder parameter 
is rather insensitive to the variation in a c . Medium to good 
quality fits can be performed for lattice sizes L > 12 for both 
models. All results for v are consistent with the best known 
value 0.7112(5) for the 3D 0(3) universality class.— Our es- 
timate for v as in Table [TV] improves the accuracy compared 
to Ref. 20 by about one order of magnitude. However, we do 
not quite reach the level presented for the bilayer models^' 4 
This could be related to the more complicated nature of the 
phase transition in planar models, where in-plane symmetries 
are broken. 

In case of the exponent p7v, good fits to Eq. d22l could be 
performed for L > 16 resulting in almost perfect agreement 
with the reference value of fi/v = 0.5188(3), which we com- 
puted from Ref. |5H Note that in case of the ladder model, 
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Figure 7: (Color online) Finite-size scaling using (a) the slope sq^ 
of the Binder parameter, (b) the staggered magnetization, and (c) the 
staggered structure factor S (n, n). These quantities are computed at 
the critical points determined in Sec, Hill 



however, the^ 2 /d.o.f increases by one order of magnitude ac- 
companied by an increase of the value of /3/v when perform- 
ing the fit at the previous estimate for a c . This indicates that 
the result of this paper indeed captures the critical point in the 
ladder model more accurately. The same observation is true 
for the exponent 77. All results for this exponent are quoted for 
lattice sizes L > 20, indicating that this quantity is harder to 
estimate. Yet, our results are still consistent or close to the ref- 
erence value. A natural check on the consistency of our results 
is a test of the (hyper)scaling relation 2/3/v = id + z - 2 + rj), 
which seems to be satisfied for nearly all cases, but it is also 
clear that rj and /3 are probably strongly correlated as they de- 
rive from almost the same quantity. 

Finally, the interested reader is referred to Ref. [23j for a 
slight extension of the current scaling analysis. In that ref- 
erence a further comparison regarding the Binder parameter 
at the critical point in different planar and bilayer Heisenberg 



Table IV: Fit results for the critical exponents v, /3/v, and 77. We 
summarize results including a variation of the critical point within 
its error bar. For the ladder model (top group of values) fit results 
and quality of fits are also given at the previous best estimate of a z . 
The bottom group are results for the plaquette model. Numbers in 
[. . .] brackets denote the^ 2 /d.o.f. For comparison relevant reference 
values for the 3D 0(3) universality class are given in the last line. 





V 


Ply" 




1.9096 -a- 


0.712(4) [1.8] 


0.516(2) [0.5] 


0.026(2) [0.2] 


1.9096 


0.711(4) [1.8] 


0.518(2) [1.1] 


0.029(5) [0.8] 


1 .9096 + (X 


0.710(4) [1.8] 


0.519(3) [2.5] 


0.032(7) [1.4] 


1.9107"* 


0.709(3) [1.7] 


0.525(8) [15.3] 


0.051(10) [12] 


1.8230 -a 


0.708(4) [0.99] 


0.515(2) [0.84] 


0.025(4) [0.15] 


1.8230 


0.706(4) [1.04] 


0.516(2) [0.40] 


0.028(3) [0.31] 


1 .8230 + (X 0.706(4) [1.10] 


0.517(2) [1.6] 


0.031(5) [0.80] 


Ref. 54 


0.7112(5) 


0.5188(3) 


0.0375(5) 



a L> 12. 
b L> 16. 
C L > 20. 

rf Previous best estimate of Ref. CO 



models is presented. 



V. SUMMARY AND CONCLUSIONS 



In conclusion, we have considered two particular geomet- 
ric arrangements of competing interactions in 2D planar quan- 
tum Heisenberg models complementing work we have started 
in Ref. 23. From detailed QMC simulations and a finite-size 
scaling study, this work provides a first high-precision value 
for the critical point in the plaquettized Heisenberg model and 
improves the value for the ladder model. In both cases, the use 
of correction terms and a combined analysis of different quan- 
tities is essential. For both models we derive the full set of crit- 
ical exponents and improve their accuracy by about one order 
of magnitude (from v = 0.71(3) to 0.711(4)) for the ladder 
model. These values are in excellent agreement with the clas- 
sical 3D 0(3) universality classi^^ As outlined above, the 
new estimates will be useful and necessary in connection with 
the recent fascinating studies on impurity based questions. In 
this regard, an extension from bilayer to planar models has yet 
to be done. 

Note added. Recently, a report by Albuquerque et al^ ap- 
peared, which also presents simulations on the plaquettized 
Heisenberg model. Since their motivation is mainly oriented 
towards showing the applicability of the contractor renormal- 
ization (CORE) method to quantum spin systems, less empha- 
sis is spent on the analysis of the critical point in detail. 



10 



Acknowledgments 

We acknowledge stimulating discussions with A. Sandvik, 
S. Wessel, and A. Muramatsu. We thank L. Bogacz for collab- 
oration in an early stage of this project. S.W. acknowledges 



support from the Studienstiftung des deutsches Volkes, the 
DFH-UFA under Contract No. CDFA-02-07, and the Leipzig 
graduate school "BuildMoNa." This work was partially per- 
formed on the JUMP computer of NIC at the Forschungszen- 
trum Jiilich under Project No. HLZ12. 



Electronic address: wenzel@itp.uni-leipzig.de 
Electronic address: janke@itp.uni-leipzig.de 
U. Schollwock, J. Richter, D. Famell, and R. Bishop (Eds.), Lec- 
ture Notes in Physics Vol. 645, (Springer, Berlin, 2004). 
S. Sachdev, Nature Phys. 4, 173 (2008). 

S. Trotzky, P. Cheinet, S. Foiling, M. Feld, U. Schnorrberger, 
A. M. Rey, A. Polkovnikov, E. A. Demler, M. D. Lukin, and 
I. Bloch, Science 319, 295 (2008). 

S. Sachdev, Quantum Phase Transitions (Cambridge University 

Press, Cambridge, 1999). 

M. Vojta, Rep. Prog. Phys. 66, 2069 (2003). 

A. W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007). 

T. Giamarchi, C. Ruegg, and O. Tchernyshyov, Nature Phys. 4, 

198 (2008). 

A. W. Sandvik and D. J. Scalapino, Phys. Rev. Lett. 72, 2777 
(1994). 

A. W. Sandvik, A. V. Chubukov, and S. Sachdev, Phys. Rev. B 51, 
16483 (1995). 

M. Troyer and S. Sachdev, Phys. Rev. Lett. 81, 5418 (1998). 

P. V. Shevchenko, A. W. Sandvik, and O. P. Sushkov, Phys. Rev. 

B 61, 3475 (2000). 

A. Collins and C. J. Hamer, Phys. Rev. B 78, 054419 (2008). 
H. G. Evertz, Adv. Phys. 52, 1 (2003). 

L. Wang, K. S. D. Beach, and A. W. Sandvik, Phys. Rev. B 73, 
014431 (2006). 

K. H. Hoglund and A. W. Sandvik, Phys. Rev. Lett. 99, 027205 

(2007) . 

S. Sachdev, C. Buragohain, and M. Vojta, Science 286, 2479 
(1999). 

F. Anfuso and S. Eggert, Phys. Rev. Lett. 96, 017204 (2006). 

O. P. Sushkov, Phys. Rev. B 62, 12135 (2000). 

M. Troyer, M. Imada, and K. Ueda, J. Phys. Soc. Jpn. 66, 2957 

(1997). 

M. Matsumoto, C. Yasuda, S. Todo, and H. Takayama, Phys. Rev. 
B 65, 014407 (2001). 

S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys. Rev. Lett. 

60, 1057 (1988). 

A. V. Chubukov, S. Sachdev, and J. Ye, Phys. Rev. B 49, 11919 
(1994). 

S. Wenzel, L. Bogacz, and W. Janke, Phys. Rev. Lett. 101, 127202 

(2008) . 

N. B. Ivanov, S. E. Kriiger, and J. Richter, Phys. Rev. B 53, 2633 
(1996). 

R. R. P. Singh, M. P. Gelfand, and D. A. Huse, Phys. Rev. Lett. 

61, 2484(1988). 

S. E. Kriiger, J. Richter, J. Schulenburg, D. J. J. Farnell, and R. F. 
Bishop, Phys. Rev. B 61, 14607 (2000). 

P. Tomczak and J. Richter, J. Phys. A: Math. Gen. 34, L461 
(2001). 

R. Darradi, J. Richter, and S. E. Kriiger, Cond. Matt. 16, 2681 
(2004). 

U. Low, Phys. Rev. B 76, 220409(R) (2007). 

N. Katoh and M. Imada, I. Phys. Soc. Jpn. 62, 3728 (1993). 

O. Nohadani, S. Wessel, and S. Haas, Phys. Rev. B 72, 024440 



(2005). 

32 O. Nohadani, S. Wessel, B. Normand, and S. Haas, Phys. Rev. B 
69, 220402(R) (2004). 

33 C. Yasuda, S. Todo, M. Matsumoto, and H. Takayama, Phys. Rev. 
B 64, 092405 (2001). 

34 T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and M. P. A. 
Fisher, I. Phys. Soc. Ipn. Suppl. 74, 1 (2005). 

35 A. Koga, S. Kumada, and N. Kawakami, J. Soc. Phy. Jpn 68, 642 
(1999). 

36 J. Sirker, A. Kliimper, and K. Hamacher, Phys. Rev. B 65, 134409 
(2002). 

37 R. R. P. Singh, Z. Weihong, C. J. Hamer, and J. Oitmaa, Phys. 
Rev. B 60, 7278 (1999). 

38 A. Lauchli, S. Wessel, and M. Sigrist, Phys. Rev. B 66, 014401 
(2002). 

39 S. Capponi, A. Lauchli, and M. Mambrini, Phys. Rev. B 70, 
104424 (2004). 

40 A. W. Sandvik and J. Kurkijiirvi, Phys. Rev. B 43, 5950 (1991). 

41 A. W. Sandvik, Phys. Rev. B 59, R14157 (1999). 

42 O. F. Syljuasen and A. W. Sandvik, Phys. Rev. E 66, 046701 

(2002) . 

43 F. Alet, S. Wessel, and M. Troyer, Phys. Rev. E 71, 036706 (2005). 

44 F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001). 

45 M. Troyer, S. Wessel, and F. Alet, Phys. Rev. Lett. 90, 120201 

(2003) . 

46 A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 63, 1 195 
(1989). 

47 M. Troyer, F. Alet, and S. Wessel, Braz. J. of Phys. 34, 377 (2004). 

48 A. W. Sandvik, Phys. Rev. B 56, 1 1678 (1997). 

49 M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, 
Phys. Rev. B 40, 546 (1989). 

50 K. S. D. Beach, L. Wang, and A. W. Sandvik, arXiv:0505194 (un- 
published) (2005). 

51 S. Wenzel, E. Bittner, W. lanke, and A. M. I. Schakel, Nucl. Phys. 
B 793, 344 (2008). 

52 K. Binder, Z. Phys. B 43, 119 (1981). 

B. Efron, The Jackknife, the Bootstrap, and other Resampling 
Plans (Society for Industrial and Applied Mathematics [SIAM], 
Philadelphia, 1982). 

54 M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi, and 
E. Vicari, Phys. Rev. B 65, 144520 (2002). 

55 C. Holm and W. Janke, Phys. Rev. B 48, 936 (1993). 

56 A. F. Albuquerque, M. Troyer, and J. Oitmaa, Phys. Rev. B 78, 
132402 (2008). 

57 For the plaquette model, there is a symmetric transition, when a is 
decreased. For the ladder model this second transition is between 
a ladder structure and Neel order. 

58 A slightly better way would be to test convergence for the largest 
lattice size and to simulate all smaller systems using the same 
temperature. 

59 For critical couplings, the method of Ref. [5(j produced reliable 
estimates which are consistent with other results in the literature. 



